earthgenome

NOAA GLOBE DEM

library(raster)

## Loading required package: sp

library(plotly)

## Loading required package: ggplot2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:graphics':
## 
##     layout

setwd("/Users/sigmamonstr/Google Drive/DOC/044-EarthGenome/sample")

file_list <- unique(substr(list.files(),1,3))

for(k in file_list){
  print(paste("Starting: ",k))
  
  #Read in header file
    hdr <- read.table(paste(k,"s.hdr",sep=""), quote="\"", stringsAsFactors=FALSE)
    hdr[,2] <-as.numeric(hdr[,2])
    col= hdr[4,2]
    row = hdr[3,2]
    mat= col*row
    xmin = hdr[11,2]
    xmax = xmin + hdr[13,2]*col
    ymin = hdr[12,2]
    ymax = ymin + hdr[14,2]*row
    print(paste("HDR: ",k))
  
  #Red in BIN file
    bil <- raster(ncols=col, nrows= row, xmn=xmin, xmx=xmax, ymn=ymin, ymx=ymax )
    bin <- readBin(paste(k,"g",sep=""), what="integer", n=mat,endian="little", signed=TRUE, size=2)
    temp <- setValues(bil, bin)
    assign(k,temp)
     print(paste("BIN: ",k))
  
}

## [1] "Starting:  a10"

## Warning: NAs introduced by coercion

## [1] "HDR:  a10"
## [1] "BIN:  a10"
## [1] "Starting:  b10"

## Warning: NAs introduced by coercion

## [1] "HDR:  b10"
## [1] "BIN:  b10"
## [1] "Starting:  e10"

## Warning: NAs introduced by coercion

## [1] "HDR:  e10"
## [1] "BIN:  e10"
## [1] "Starting:  f10"

## Warning: NAs introduced by coercion

## [1] "HDR:  f10"
## [1] "BIN:  f10"

  #Merge US rasters
    north_america <- merge(e10,f10)
  
  #Crop to continental 
    e <- extent(-130, -60, 70, 100)
    north_america <- crop(north_america, e)  

  #Average every 10 rasters
    r <- aggregate(north_america, fact=10, fun=mean)
    r <- as.matrix(r)
    a<-r

  #Center and standardize raster values
    a <- (a)/median(a)
    a[a < 0]<-NA

  #Flip order of columns  
    a<-a[,seq(dim(a)[2],1,-1)]

  #Parameter for view of 3d surface
    scene = list(
        camera = list(
            eye = list(x = 0, y = 1, z = 1)
        ),
        aspectratio =  list(x = 1, y = 1, z = 0.3)
        )
    
        
    
      plot_ly(z = a, type = "surface", 
              colors = terrain.colors(10)) %>% 
        layout(scene=scene)